Probability of metastable configurations in spherical three-dimensional Yukawa 

crystals 
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Recently the occurrence probabilities of ground- and metastable states of three-dimensional 
Yukawa clusters with 27 and 31 particles have been analyzed in dusty plasma experiments [Block 
et al.. Physics of Plasmas 15, 040701 (2008)]. There it was found that, in many cases, the ground 
state appeared substantially less frequently than excited states. Here we analyze this question the- 
oretically by means of molecular dynamics (MD) and Monte Carlo simulations and an analytical 
method based on the canonical partition function. We confirm that metastable states can occur 
with a significantly higher probability than the ground state. The results strongly depend on the 
screening parameter of the Yukawa interaction and the damping coefficient used in the MD simula- 
tions. The analytical method allows one to gain insight into the mechanisms being responsible for 
the occurrence probabilities of metastable states in strongly correlated finite systems. 



I. INTRODUCTION 

Finite strongly coupled systems of charged particles in 
external traps are of high interest in many fields. Ex- 
amples include ion crystals [l], Q , quantum dots |3|] and 
dusty plasma crystals [1, . Dusty plasmas allow for an 
easy realization of strong coupling in laboratory experi- 
ments. They typically consist of /im sized particles in an 
rf discharge. Due to their high mass their motion occurs 
on a macroscopic timescale which makes them an ideal 
system for studying dynamical properties in the strong 
coupling limit. In the case of an isotropic parabolic con- 
finement and (screened) Coulomb interaction the ground 
states are nested spherical shells (3D) or concentric rings 
(2D). 

For classical systems the ground states are found by 
minimizing the potential energy U with respect to all 
particle positions. This can be a difficult task since in 
general U has many minima which may be energetically 
very close to each other, particularly in 3D. To find the 
lowest energy configuration one has to avoid trapping 
in a metastable state, which can be a serious problem 
for numerical computations. A detailed analysis of the 
ground states of 3D Coulomb clusters was presented in 
E] . The ground states of small spherical Yukawa clusters 
for a wide range of the screening parameter can be found 
in (sj. Besides the ground state also metastable states 
were found in the simulations d, 0) Furthermore a 
fine structure was observed, i.e. states with the same 
number of particles on each shell but with a different 
arrangement on the same shell @ . 

Coulomb or Yukawa balls have been produced in dusty 
plasma experiments They are well explained by a 
simple model of harmonically confined particles inter- 
acting by a Yukawa potential for TV = 100 . . . 500 [lo| . 
Recently metastable states of Yukawa balls have been in- 



vestigated in [llj for small particle numbers TV = 27 and 
31. It was found that often metastable states oc- 
curred with a higher probability than the ground state. 
This was confirmed by MD simulations but no theoreti- 
cal explanation was given. This is the goal of the present 
paper. We apply Monte Carlo simulations (MC) as well 
as extensive molecular dynamics (MD) simulations with 
a broader parameter range than before, confirming the 
main results of [ll| . For a theoretical explanation we ap- 
ply an analytical method based on the classical canonical 
partition function [T2| . 

This paper is organized as follows. In Sec.[n]we present 
the Hamiltonian and explain our simulation methods. 
Results of the MD simulations are given in Sec. IIIIl In 
Sec. llVl we introduce an analytical method for the prob- 
abilities of stationary states in thermodynamic equilib- 
rium. The results are compared to MC simulations. Sec- 
tion |V] compares the theoretical results with the exper- 
iments. The last section summarizes our findings and 
discusses the applicability range of our models. 



II. MODEL AND SIMULATION IDEA 

A. Hamiltonian 

The system of N identical particles with charge Q and 
mass m in an isotropic, parabolic confinement 



(r = \r\) is described by the Hamiltonian 
TV , 2 ~i ^ 

i = l ^ ^ i>j 



(1) 



(2) 
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The interaction is assumed to be a shielded Coulomb- 
potential of the form 



V{t) 



(3) 



where the range of the potential is controlled by the 
screening parameter n. Despite its simplicity this model 
is of relevance for many systems, such as colloids, and has 
proven to accurately describe the spherical dust crystals 
(Yukawa balls) observed in experiments [l^]- In dusty 
plasmas k is given by the inverse Debye screening length. 
In the following we will treat it as a free parameter and 
focus on the general behavior of the model ^ . Using, as 
a particular example, typical dusty plasma parameters 
will allow us to make comparisons with the experimental 
observations of Ref . [ll| . 

Results will be given in units of the distance rg = 
i^Q"^ Imio^^l'^ and the corresponding Coulomb energy 
^0 = /fa- Frequencies and forces are given in units of 
loq and muj^rQ, respectively. 

The ground (metastable) states are the global (local) 
minima of the potential energy J7, 



N 



N 



i>3 



In both cases the total force on all particles vanishes and 
the system is in a stable configuration, i.e. stable against 
small perturbations. 



B. Monte Carlo (MC) 

The MC simulations use the standard Metropolis al- 
gorithm [l^ with the Hamiltonian but without the 
kinetic energy part. Starting from the classical ground 
state at T = the system is given a finite temperature. 
For a fixed temperature we performed 10^ MC steps and 
determined the configuration every IQHh step. The tem- 
perature is then increased and the same procedure re- 
peated. Ergodicity of the procedure was checked by us- 
ing different initial configurations. Following this method 
we calculate the probability as a function of T from the 
number of occurrences of the different states. 



C. Molecular dynamics (MD) 

In the MD simulations we follow a different approach. 
Here we solve the equations of motion for particles in 
a parabolic trap interacting through the Yukawa poten- 
tial ([3]) but include an additional damping term to sim- 
ulate the annealing process the way it occurs in the ex- 
periment, as explained in llj. This is different from the 
MC simulations where the particles are in contact with 
a heat bath and maintain a constant temperature. This 
also differs from the MD simulations in which were 



also performed at finite temperature. Here, we perform 
substantially larger simulations and systematically scan 
a broad parameter range. For the i-th particle the equa- 
tion of motion we solve is 



(5) 



where v is the collision frequency which will be given in 
units of uq. In dusty plasmas friction is mainly due to 
the neutral gas. 

The simulation is initialized with random particle po- 
sitions and velocities in a square box. To stop the simula- 
tion and determine the configuration we use two similar, 
but not equivalent conditions: 

(A) The particles' mean kinetic energy drops below a 
value (EJ^^) of typically 10"^ - 10"^. 

(B) The force on each particle due to the confinement 
and the other particles decreases below 10"''. 

It is tempting to define (A) as a proper condition but we 
will show that (B) has to be used, although they look 
equivalent at first glance. The difference lies in the defi- 
nition of a stable configuration. If the particles lose their 
initial kinetic energy before they have reached a local 
minimum the simulation could be stopped before the par- 
ticle motion has effectively ended. This problem can be 
circumvented by condition (B) which makes direct use of 
the definition of a stable state, namely that the force on 
each particle due to U vanishes. 

The screening parameter, the friction coefficient as well 
as the lower limit for the mean kinetic energy are varied. 
For each parameter setting the simulation is repeated 
3000 — 5000 times to obtain accurate statistics. We con- 
sider systems with 31 and 27 particles as was done in the 
experiment. As another example we used a cluster with 
40 particles because here the ground state shell configura- 
tion abruptly changes from (34,6) to (32,8) at k = 0.415 
as the screening parameter is increased - without the con- 
figuration (33,7) ever being the ground state :8i|. This 
gives rise to the question of how often this configuration 
can actually occur in experiments. 



III. MD SIMULATION RESULTS 

In this section we present the results of our first- 
principle MD simulations. The main parameters deter- 
mining the occurrence frequencies of different metastable 
states for a given N are the screening parameter k and 
the friction i/. We therefore discuss the dependence on 
K and v in detail. As an example of particular interest 
we will consider the parameter values of dusty plasma 
experiments which are in the range of k « 0.4 . . . 1.0 
This case will be dealt with in Section IVl 

We first discuss the effect of the damping rate on the 
occurrence probabilities. It will turn out that with a 
properly chosen rate we can produce very general results 
for different screening lengths which do not depend on 
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TABLE I: Energy diflterence between metastable states and 
the ground state (the ground state and its energy is given by 
italic numbers) as seen in Fig. [T] States with the same shell 
configuration but different energy differ only by the arrange- 
ment of the particles on the same shell (fine structure). 
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config. 




{27,4) 


0.000479 


(26,5) 


0.000006 


(27,4) 


0.000499 


(26,5) 


0.000009 


(27,4) 


0.000530 


(26,5) 


0.000291 


(26,5) 


0.000656 


(25,6) 


0.000372 


(26,5) 


0.000669 


(25,6) 



FIG. 1: (Color online) Stationary states observed in the MD- 
simulations for iV = 31, k = 1.4 and (^Sn") = 10"^ The 
runs are sorted by the energy or the stationary state, see 
also table m For slow cooling (black bars, v = 0.05) one can 
clearly see distinct states which correspond to the horizon- 
tal lines. The length of the bold lines is proportional to the 
occurrence probabifities. In the case of strong friction (red, 
dashed line, v = 5.3) the particles often lose their kinetic en- 
ergy before they can settle into the equilibrium positions and 
the fine structure (different states with same shell configura- 
tion) cannot be resolved. 

the exact chosen damping coefficient and hold for any 
rate in the overdampcd hmit. The effect of screening 
will then be examined in the following section. 

A. Effect of friction 

A typical simulation result is shown in Fig. [TJ For 
slow cooling (y = 0.05) the particles are not hindered 
by friction and can move according to the interparticle- 
and confinement forces. They continuously lose kinetic 
energy until they are trapped in a local minimum of the 
potential energy U . Here they are further being damped 
until the simulation is stopped. It is interesting to see 
that there exist more metastable states than different 
shell occupations, as was first observed in see also f^. 
Details are given in Table [J 

In the case of strong damping [y = 5.3) the situation 
is different. Here the particles are readily slowed down 
after the initialization process in the box. Their mo- 
tion is strongly affected by friction and interrupted even 
before they may be trapped in a local minimum. If con- 
dition (A) is used to stop the simulation it is not clear if 
the particles are in a stable state. The reason is that due 
to the rapid damping they can be sufficiently slowed even 
though they are not in a potential minimum but on a de- 
scending path and would reach the stable configuration 
at a later time. 

Fig. [2] shows the influence of friction on the occurrence 
probabilities in more detail. For fixed screening the prob- 
ability of finding the ground state configuration increases 
when the friction coefficient is decreased. Here the par- 



ticles are cooled down more slowly and it is more likely 
that they reach the system's true ground state. During 
the cooling process they still have a sufficiently high ki- 
netic energy and time to escape from a local minimum 
until the force on each particle vanishes. In the case of 
strong friction the particles can fall into a nearby min- 
imum and leaving it becomes more difficult due to the 
rapid loss of kinetic energy. The typical simulation time 
until the forces are small enough is longer than for inter- 
mediate friction strength. Once cooled down the parti- 
cles are pushed along the gradient of the potential energy 
surface until they reach a stable state. Thus the results 
can depend on how far the system's temperature is de- 
creased. One can see that for > 2, i.e. in the over- 
damped regime, the probabilities have practically satu- 
rated. For fast cooling, i.e. large friction, metastable 
states can occur with a comparable or even higher prob- 
ability than the ground state. 

The = 40 cluster shows a qualitatively different be- 
havior compared to the A^ = 27, 31 clusters. In the case 
of K = 1.0 the lines corresponding to different configura- 
tions do not intersect and the ground state is the most 
probable state regardless of the damping coefficient. In 
contrast, in the Coulomb limit, k = 0, the most probable 
state is a always a metastable state, except for very small 
friction, v < 0.01. 

Dusty plasma experiments are performed in the over- 
damped regime, i.e. here v is of the order of 3 — 6 jTl| . 
Since in this limit the probabilities depend only very 
weakly on the damping rate the results presented in the 
next section for v = 3.2 should hold for any such damp- 
ing coefficient. Even though this was shown only for a 
few examples we believe that this also holds for other 
particle numbers and screening lengths. 



B. Effect of screening 

The screening dependence of the ground state shell 
configurations of spherical Yukawa clusters in the absence 
of damping has been analyzed in Ref. p^ . The general 
trend is that increased screening favors ground state con- 
figurations with more particles on the inner shell(s). A 
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a) N = 27, K = 0.6, groundstate: (24,3) 




Qg_ b)N = 31, K = 0.8, groundstate: (27,4) - 
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FIG. 2: (Color online) Effect of friction on the occurrence 
probabilities obtained with condition (B) for three different 
numbers of particles. In a) and b) horizontal solid and dashed 
lines indicate experimental mean and standard deviation, re- 
spectively [13]. For N = 27 the experimental values for the 
clusters (23,4) and (24,3) are the same. In c) solid lines in- 
dicate Yukawa interaction with k = 1.0 [ground state (32,8)] 
whereas dashed lines show results for Coulomb interaction 
[ground state (34,6)]. In all cases slow cooling favors the 
ground states over metastable states. 
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FIG. 3: (Color online) Effect of screening for f = 3.2. Solid 
lines show the results obtained with condition (B) while 
dotted and dashed lines indicate use of condition (A) with 
(^'Sri") ~ 10~*, 10"^, respectively. Arrows show the ground 
state configuration to the left or right from the vertical line. 
Where available horizontal solid and dashed lines indicate ex- 
perimental mean and standard deviation [ll|] . For the N — 27 
cluster the experimental values for the configurations (23,4) 
and (24,3) are the same. 



systematic analysis in a large range of particle numbers 
and screening parameters 8] confirms this trend. Here we 
extend this analysis to spherical crystals in the presence 
of damping and also consider the screening dependence 
of the occurrence probability of metastable states. 

For a fixed friction coefficient in the overdamped limit 
the effect of screening is shown in Fig. [31 The different 
ground state configurations are indicated by the numbers 



with arrows in the figures. As in the undamped case, at 
some finite value of k, a configuration with an additional 
particle on the inner shell becomes the ground state. 
Consider now the probability to observe the ground and 
metastable states. For weak screening the ground states 
(27,4) and (24,3) are the most probable states in the cases 
iV = 31 and N = 27, respectively. At the same time in 
both cases, the probability of the configuration with one 



5 



more particle on the inner shell grows with k, until it 
eventually becomes even more probable than the ground 
state. Note that this occurs much earlier (at a signifi- 
cantly smaller value of k) than the ground state change. 
For = 31 this trend is observed twice: the proba- 
bility of the configuration (26,5) first increases with k 
and reaches a maximum around k w 1. For k > 2 this 
configuration becomes less probable than the configura- 
tion (25,6), i.e. again a configuration with an additional 
particle on the inner shell becomes more probable with 
increased screening. 

Different behavior is observed for the = 40 cluster 
where the ground state for weak screening (34,6) is never 
the most probable state. For large screening, k > 0.6, the 
new ground state (32,8) has the highest probability, but 
this happens only substantially later (for larger k) after 
this state has become the energetically lowest one. This 
is due to the existence of a third state (33,7) which has 
the highest probability for /t < 0.6 although it is never 
the energetically lowest one. 

Summarizing the above observations we confirm that 
in spherical Yukawa clusters the ground state is not nec- 
essarily the most probable state. Often, a metastable 
state with more particles on the inner shell is observed 
substantially (in some cases up to five times) more fre- 
quently. Further, increased screening tends to favor 
states with more particles on the inner shell. We will 
give an explanation for this behavior in the next section 
by using an analytical model for the partition function. 

Before doing this we comment on some technical de- 
tails which are important in the present MD simula- 
tions. For certain intervals of the screening parameter 
the results for the probabilities depend on how far the 
system is cooled down. Here one state (generally the 
ground state) is favored over another the smaller {E'^"^) 
is chosen. This also means increasing the mean simula- 
tion time. As discussed before the particles are heavily 
damped and lose their initial kinetic energy on a short 
timescale. Their motion is then determined by the shape 
of the energy surface. Using condition (B) to terminate 
the simulation we obtain converged results where the par- 
ticles have reached a local minimum. Thus if the simula- 
tion would be continued the configuration would remain 
the same. 



IV. ANALYTICAL THEORY OF STATIONARY 
STATE PROBABILITIES 



the potential energy of a given state is expanded around 
a local minimum with energy E^, where s denotes the 
ground- or a metastable state. It can be written as 



1 



N 



2 dri^adrj^fs 

i,J = l a, 13=1 ' -"^ 



Svi.aSr 



where r"" = (rj'', . . . , r^) denotes the 3Af-dimensional 
vector of the particles' equilibrium positions and 5ri a ~ 
i'i,a ~'"?a the displacement vector. Transforming to nor- 
mal coordinates ^s,i this turns into a sum of decoupled 
harmonic oscillators 



1 



1=1 



(6) 



with eigenfrequencies LOs,i, which are the square roots of 
the eigenvalues of the Hessian 



i,a,j,f3 



'3,13 



The expansion ([6]) includes the particles' three center of 
mass oscillations in the trap with w = 1 (in units of 
loq). Furthermore we assume that the vibrational and 
the three rotational modes of the whole system (w = 0) 
are decoupled, the latter are, therefore, eliminated from 
the sum In the principal axes frame the rotational 
kinetic energy can then be expressed as 



E 



2Ia_i 



with angular momenta Ls^i and constant principal mo- 
ments of inertia Ig^i- In this approximation the full en- 
ergy of the state s is, to second order in the displace- 
ments, 



E., 




3 r2 



i=l 



2L 



(7) 



The first term in parentheses denotes the vibrational ki- 
netic energy TJ*''. 

The harmonic approximation is only applicable for low 
temperatures (or strong coupling) when the particles os- 
cillate around the equilibrium positions with a small am- 
plitude. 



A. Harmonic approximation 



Partition function 



The analytical approach to calculating the occurrence 
probabilities is based on the classical canonical partition 
function Z{T,ujo, N). Instead of the dependence on vol- 
ume (or density) as in a homogeneous system, here ther- 
modynamic quantities depend on the confining strength 
uq. The partition function can be evaluated analytically 
in the harmonic approximation, see e.g. Ref. }12\. Here 



The general form of the classical canonical partition 
function is 



(271^)3^ J_ 



(8) 



Here it is written for a general Hamiltonian H^{pi,qi) 
with iN degrees of freedom, generalized coordinates qi 
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and conjugate momenta pi. Since in our case the en- 
ergy contributions are independent it can be factorized 
according to 



Ze TleZ t" ' Z', 



ib ^^ot 



with the internal partition function 



and the degeneracy factor Ug calculated as 

m 



(9) 



(10) 



(11) 



where L is the number of shells and Nf the occupation 

number of shell i with X^i^Li -^i = The degeneracy 
factor Hs denotes the number of possibilities to form 
a configuration with shell occupation (A^i, N2, . . . , Nl) 
from distinguishable particles. 

ZJ*'' is the partition function for / independent har- 
monic oscillators while ZJ"* is related to the rotational 
degrees of freedom. The results for our specific case with 
the energy given by Eq. ([7]) can be found in and read 



Mis 



3/2 



(12a) 
(12b) 



The expressions include the mean geometric eigenfre- 
quency VLs = (I\i=i^s,iY^ ^ and the mean moment of 
inertia Is = {Is,iIsaIs,?.Y^^ ■ 

To obtain the total partition function Z{T,u!o^ N) the 
contributions of all M (metastable) states are summed 
up, i.e. 

M 



Z = ^ n^Za 



C. Probability of stationary states 



Collecting the results of subsection llV B[ the stationary 
state probabilities are given by 



Ps 



Us Zs 



TlsZs 



(13) 



For our clusters of interest with 27 — 40 particles the mo- 
ments of inertia for different states are equal to a good 
approximation (cf. Table ITTl for TV = 27) and can be can- 
celed. Similar behavior is observed for — 31, 40. For 
low particle numbers, N < 10, they should be included, 
since here a slight change of the configuration can alter 
the moment of inertia by a significant amount, but this 
is not of importance for the present analysis. 



TABLE II: Mean shell radii _Ri, R2 of first and second shell 
for states observed in the MD simulations for N — 27 and 
K = 0.6. The relative statistical weight Qs = (Is/Ii)^^^ caused 
by different moments of inertia can be neglected in the com- 
putation of the probabihties since Qa ~ 1 for all states. 



state s 


configuration 


R2 


Ri 




1 


(24,3) 


1.6175 


0.5977 


1 


2 


(23,4) 


1.6413 


0.6963 


1.0009 


3 


(23,4) 


1.6413 


0.6957 


1.0009 


4 


(25,2) 


1.5935 


0.4542 


1.0004 


5 


(25,2) 


1.5934 


0.4543 


1.0004 



Using Eqs. we obtain from Eq. ^3]) 
nse-P^"n-f 



(14) 



To avoid computation of the full partition function [de- 
nominator of Eq. (|14p ] it is advantageous to compute 
probability ratios of two states s and s' 



Ps' 



Us 
Hs' 



^s 

Qs^ 

^s 



(15) 



Thus the probability ratio of two states depends on three 
factors: their energy difference — , the ratio of de- 
generacy factors Us/ Us' and the ratio of mean eigenfre- 
quencies il^z/fJ^.. 

The Boltzmann factor e^^(^?~^a') gives preference to 
states with a low energy. For low temperatures it will 
be the most dominant factor but it becomes less impor- 
tant for higher temperatures when /c^T 3> E^, — and 

g-/3(B!?-i5«,) _ ^ 

According to Eq. pT|) the degeneracy factor assigns 
a large statistical weight to states with more particles 
on inner shells. As an example, for N = 27, we ob- 
tain ?^(25,2)/'^'(23,4) = Hrfl = 1/50. One can see that the 
configuration with only 2 particles on the inner shell is 
suppressed due to a lower degeneracy factor contrary to 
the states with an inner shell consisting of 4 particles, see 
also Table IIIIl The reason is that there exist more com- 
binatorial possibilities to construct configurations when 
the difference between the single shell occupation num- 
bers is small. For = 31 (Table |V| this ratio can be 
even larger. This shows that (even for low temperatures) 
this factor can strongly influence the occurrence proba- 
bilities. 

In the MD simulations we observe several states with 
the same shell configuration but different energies. Their 
energy difference can be as large as between states with 
different configurations (cf. Table. |l|. In Eq. (fT3|) all 
states with the same shell configuration are added with 
the same degeneracy factor. 
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Let us now consider the effect of tfie mean eigenfre- 
quency, i.e. the effect of the local curvature of the poten- 
tial energy surface. Written out explicitly, using Eq. (fTS]) . 
this factor reads 



nf=i ^s,^ ' 



(16) 



i.e. it is the inverse ratio of the products of the eigenfre- 
quencies. The main contribution here usually arises from 
the lowest eigenfrequencies. This can be seen in Fig. |3] 
showing the spectrum for the states of the cluster with 
N = 31, K — 0.8. State #7 has two very low eigenfre- 
quencies [cf. Fig. |4l red arrow] which strongly increase 
its statistical weight (see also Table ITVT) . 
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FIG. 4: Spectrum of the eigenfrequencies for the 9 states 
shown in Table ITVl The top figure shows the lowest modes in 
more detail. 

For two states with the same shell configuration we 
have Us = Ugi , and the probability ratio is only deter- 
mined by their energy difference and eigenfrequencies. 
Even though a state has a higher energy it can have a 
higher probability provided it has a lower mean eigenfre- 
quency. Fig. [5]shows the effect for N = 27, for the states 
listed in Table Hill The physical explanation of the eigen- 
frequency factor is very simple: states with low eigenfre- 
quencies have a broad (flat) potential energy minimum 
and thus a larger phase space volume of attraction for 
the trajectories of N particles. Thus initially randomly 
distributed particles will have a higher probability to set- 
tle in a minimum with small compared to another 
minimum (when the energies and degeneracy factors are 
similar). 

Because the harmonic approximation only describes 
a minimum's local neighborhood we mention that this 
could overestimate the weight of states with broad min- 
ima and low escape paths [12l |. which are not taken into 



TABLE III: Energy difference between metastable states and 
the ground state (ground state energy given in italic numbers) 
that were used to compute the partition function for N — 27 
and K = 0.6. Also shown is the relative statistical weight — 
Us/ni and the statistical weight due to the eigenfrequencies 
Ws = (fii/fis)^ compared to the ground state. 



state s 


configuration 


AEs/N 




Ws 


1 


(24,3) 


4.732856(4) 


1 


1 


2 


(23,4) 


0.001622(1) 


6 


0.24 


3 


(23,4) 


0.001870(5) 


6 


0.67 


4 


(25,2) 


0.004993(0) 


3/25 


14 
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(25,2) 


0.004997(3) 


3/25 


3.3 
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FIG. 5: Probability of the two metastable states with conf. 
(23,4) compared to the ground state (24,3) for the Yukawa ball 
with = 27. The inset shows the ratio of the probabilities for 
states 2 and 3 from Table UlIl at low temperatures. Although 
state 3 has the same configuration and a higher energy the 
probability of finding state 3 is higher for T > 0.007 due to 
the effect of the eigenfrequencies. 



account in this approximation. This could be improved 
by changing the limits for the position integration in 
Eq. ([5]) according to the potential barrier height and the 
temperature. This was done for 2D clusters in [IJ] but 
requires knowledge of the barrier heights. This is not es- 
sential for the present analysis. Finally, we note that the 
value of Ws is sensitive to numerical errors in the compu- 
tation of the eigenvalues of the Hessians since the mean 
eigenfrequency is a product of 3A'^ — 3 single values. In 
the present results we estimate the error not to exceed 
5 % which is sufficient for our analysis. 



D. Analytical results and comparison with Monte 
Carlo simulations 

Let us now come to the results of the analytical model 
and compare them to Monte-Carlo simulations which 
were explained above in Section III Bl The MC results 
have first principle character, in particular, they are not 
restricted to the harmonic approximation and fully in- 



TABLE IV; Same as Table Im] for A/' = 31 and k = O.i 



ite s 


configuration 


AE°/N 


fls 


Ws 


1 


(27,4) 


4.397858(8) 


1 


1 


2 


(27,4) 


0.000008(7) 


1 


0.82 


3 


(27,4) 


0.000035(8) 


1 


1.7 


4 


(26,5) 


0.001810(1) 


27/5 


0.84 


5 


(26,5) 


0.001850(9) 


27/5 


1.4 


6 


(26,5) 


0.002000(0) 


27/5 


5.3 


7 


(26,5) 


0.002091(6) 


27/5 


9.7 


8 


(25,6) 


0.003583(7) 


117/5 


1.4 


9 


(25,6) 


0.003586(7) 


117/5 


1.1 



TABLE V: Same as Table|TlT]for A = 40 and k = (Coulomb 
interaction) . 



ite s 


configuration 


AE°/N 




Ws 


1 


(34,6) 


12.150162(9) 


1 


1 


2 


(33,7) 


0.001143(4) 


34/7 


2.3 


3 


(33,7) 


0.001190(3) 


34/7 


2.8 


4 


(33,7) 


0.001236(9) 


34/7 


8.3 


5 


(32,8) 


0.001862(8) 


561/28 


13 


6 


(32,8) 


0.001863(1) 


561/28 


3.5 


7 


(32,8) 


0.003482(4) 


561/28 


6.7 


8 


(35,5) 


0.004201(7) 


6/35 


5.2 


9 


(35,5) 


0.004392(7) 


6/35 


32 



elude all anharmonie corrections. For N = 27 we addi- 
tionally verified the MC results by a Langevin dynamics 
simulation using the SLO algorithm of [1^1 ■ Here the 
probabilities were obtained in an equilibrium calculation 
with a simulation time t = 10^ luq^ by determining the 
configurations at fixed time intervals. 

Results for three representative examples are shown in 
Fig. M We chose N ^27, k = 0.6 and TV = 31, k = 0.8 
since these will turn out to be close to the situation in 
the dusty plasma experiments, see. Sec. |Vl As a third 
example we present data for A'^ = 40 with Coulomb inter- 
action. The input parameters of the analytical model, i.e. 
details on the (metastable) states are summarized in Ta- 
bles IIIHTVl In Fig. [5] we plot the occurrence probabilities 
as a function of temperature. This allows us to specifi- 
cally study the effect of tlie deptli of tlic potential energy 
minimum Es. The latter effect should be dominant at 
low temperature, leading to a relatively high probabil- 
ity of the ground state. In contrast, this effect should 
become less important at high temperature where the 
degeneracy factors and the eigenfrequency ratio should 
play a decisive role for the probabilities. This general 
trend is indeed observed in all three cases. 

For N — 27, top part of Fig. [6l the effects of the de- 
generacy factor and the mean eigenfrequencies act in op- 
posite directions. While the state with 4 particles on 
the inner shell gains statistical weight by having a high 
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FIG. 6: (Color online) Analytical theory compared to MC- 
results. The solid lines show the probabilities as obtained 
from Eq. (I13|l . The dashed lines neglect the statistical weight 
factor caused by the eigenfrequencies, i.e. here fis = 1 for 
all states. For N = 27 the dashed/dotted lines indicate the 
results of Langevin dynamics simulations. Analytical results 
for the configuration (22,5) are not available. 



degeneracy, this effect is almost compensated by nar- 
row minima and, consequently, a low Ws, cf. Tab. IIIII 
Therefore, this state achieves comparable probability 
with the ground state (27,4) only at high temperature, 
T > 0.03 (in the MC simulation this is observed only for 
T > 0.045). For the configuration (25,2) the opposite is 
true. Here, the degeneracy is low and the minima broad, 
but due to its high energy this configuration has a nonva- 
nishing probability only for high temperatures, T > 0.03. 
We did not find a stable state with configuration (22,5). 
The situation for = 31, central part of Fig. [51 is 
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different. Here all metastable states have a higher de- 
generacy factor than the ground state configuration. In 
addition all states further gain statistical weight because 
of broad minima, except for state s = 4, cf. Tab. IIVI 
Thus one should expect that metastable states have a 
high probability even at low temperatures. This is indeed 
observed in the model and the MC simulations already 
below T = 0.02. 

In the third case, N = 40, bottom part of Fig. [6l we 
generally see the same trend. The metastable state (33,7) 
has a high degeneracy and frequency factor, cf. Tab. IIVI 
and thus it becomes more probable than the ground state 
already for T > 0.01 (0.015 in the MC simulations). 

Let us now compare the analytical and MC results 
more in detail. Good agreement is found for = 31 up to 
T ~ 0.02, cf. full lines and symbols. For TV = 27 wc find 
good agreement between MC and the analytical theory 
for T < 0.012 but only if the effect of the eigenfrcqucncies 
is neglected, cf. dashed lines. With eigenfrequencies in- 
cluded the theory shows deviations for low temperatures 
but better agreement for higher temperatures. For the 
cluster with 40 particles we observe moderate agreement 
for the configurations (34,6) and (33,7) up to T = 0.015 
whereas the deviations from MC for the remaining two 
configurations are rather large. This overall agreement is 
quite satisfactory keeping in mind that the melting tem- 
perature of these clusters is typically below T = 0.015 

til,[i3,[a. 

The reason for these discrepancies are due to the 
limitations of our simple harmonic model [the good 
agreement between the completely independent MC and 
Langevin MD results for N — 27, cf. top part of Fig. [51 
confirms the reliability of the simulations] . Since the dis- 
crepancies are growing with temperature, the main rea- 
son is probably the neglect of anharmonic effects. In 
some cases, when the barriers of the potential energy sur- 
face are low, these effects might already occur at low tem- 
peratures. Changing the limits of allowed particle motion 
in the integration of Eq. ^ may help to reduce the devi- 
ations. A further reason for deviations from MC results 
could be an insufficient number of stationary states being 
taken into account. It is not clear if all stationary states 
have been found (they were pre-computed with MD sim- 
ulations) and used in the partition fimction. To ensure 
a high probability we performed more than 10* indepen- 
dent runs. For example, for the cluster with 40 particles 
we observe 9 states, but it was difficult to identify the 
states with 5 particles on the inner shell because they 
were found only a few times and were energetically close. 
The larger number of states given in [9| also suggests 
that we missed a few. Nevertheless, the effect originat- 
ing from these states should give only a small statistical 
contribution to the probabilities. 



V. COMPARISON WITH DUSTY PLASMA 
EXPERIMENTS 

To compare with experiments on metastable states in 
spherical dusty plasma crystals (Yukawa balls) we first 
need to establish the relation of our system of units 
to the experimental parameters. We use the tempera- 
ture unit keTo = Eq = {aQ^ [in SI units = 
(aQ'^ l'yiir^(^)^l'^\ which depends on the trap parame- 
ter a. = towq and the dust charge. Since the charge 
is not known very accurately the errors could be rather 
large. With Z = 2000 e and a = 5.2 x 10""kgs-2 
given in [TT'| , room temperature (300 K) corresponds to 

oorn ~ 0.0015. Also, the experimental screening pa- 
rameter is known only approximately. From previous 
comparisons with simulations ^] it is expected to be 
in the range of 0.5 < k < 1. Reference [ll| reported 
measurements on the probability of metastable states for 
two clusters with N — 27 and = 31 which we now use 
for comparison with the MD and MC simulations and 
the analytical model. 



A. MD results vs. experiment 

We start with the molecular dynamics simulations 
since they model a situation which is closest to the ex- 
periment. In contrast to the experiment which is per- 
formed at room temperature, our simulations correspond 
to a Langevin dynamics simulation at T = (the sys- 
tem is cooled to almost zero kinetic energy). We have 
verified the influence of the final temperature by per- 
forming additional Langevin simulations for the cluster 
with 27 particles and n = 0.6 with temperatures up to 
T — 0.0035 (Fig. [7]) which is more than twice the exper- 
imental temperature. Apart from a finite temperature 
the simulations were done in the same way as explained 
in Section [ni but with a predefined simulation time. For 
high temperatures one has to pay attention to the time 
after which the configuration is determined since then 
transitions between states can easily occur. This can 
be seen in Fig. [5] where for T > 0.01 metastable states 
have a nonvanishing probability. In our Langevin sim- 
ulations we used a simulation time of t^nd = 400 w,^^, 
which corresponds to tend ~ 10 s for a dust particle mass 
of TO = 3.3 X 10~^*kg. We find no systematic deviation 
from the results at zero temperature. The slight devia- 
tions for the configurations (23,4) and (24,3) are probably 
due to the insufficiently long simulation time with the 
same explanation as given at the end of Section IIIIBI 
We thus conclude that for the present analysis an MD 
simulation without fluctuations and cooling towards zero 
temperature is adequate. 

Our data for comparison with the experimental results 
are shown in Figs. [2] and [3l The friction parameter in 
the experiments is expected to be in the range ly — 3 . . .6 
[m . This means the system is overdamped and any value 
above v = 2 will not change the results signiflcantly, cf 
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FIG. 7: Langevin dynamics simulation for = 27, k — 0.6 
and u = 3.2. Horizontal lines indicate results of Section FlII Bl 
Fig. El 

Fig. H So in Fig. [3] we used a value of 3.2. The MD 
simulations agree well with the experiment in the case 
of screening parameters in the range 0.6 < k < 0.8 (for 
N = 31) and 0.4 < k < 0.6 (A^ = 27), for details cf. Ta- 
ble |Vll The lower screening parameter in the latter case 
is a consequence of the lower plasma density in the exper- 
iment, compared to the conditions under which the clus- 
ter with 31 particles was produced. This was also found 
in the MD simulations performed in [ll[. The present 
simulations, being much more extensive, confirm these 
results. We may conclude that this comparison allows to 
determine the screening parameter in the experiment. 

B. Analytical and MC results vs. experiment 

A comparison of the analytical model and the MC 
simulations with the experiment is disappointing. From 
Fig. [5] it is evident that at room temperature the ground 
states have always a probability of almost 100 % which 
is in striking contrast to the experiment and the MD re- 
sults. This is not surprising since the dust comprises 
a dissipative system and the clusters are created un- 
der nonequilibrium conditions. In contrast, both Monte 
Carlo and the model are based on the canonical partition 
function and assume thermodynamic equilibrium. Thus, 
at first sight, there seems to be no way to explain the ex- 
periment with our analytical model or with Monte Carlo 
methods. However, this is not true. As we will show be- 
low, there is a way to apply equilibrium methods to the 
problem of metastable states. 



C. Time scales of the cluster dynamics 

Let us have a closer look at the nonequilibrium dy- 
namics of the cluster during the cooling process. It is 
particularly interesting to analyze on what time scales 



the different relaxation processes occur. In a weakly cou- 
pled plasma there are three main time scales, e.g. [l9.l20|: 
first, the buildup of binary correlations which occurs for 
times shorter than the correlation time Tcor- Second, the 
relaxation of the velocity distribution towards local equi- 
librium due to collisions, for Tcor l£ t < trei (kinetic 
phase) and third, hydrodynamic relaxation, trei < t < 
thyd- This behavior has so far not been analyzed for the 
strongly correlated Yukawa clusters. 

To get first insight, the quantities of central interest 
are the kinetic energy and the velocity distribution func- 
tion f{v,t) of the cluster particles. These quantities are 
easily computed in our nonequilibrium MD simulations 
of the coohng process, as explained in Section [1111 To ob- 
tain the velocity distribution we performed 420 runs with 
different randomly chosen initial conditions and collected 
the data for each time step. The results for the kinetic 
energy evolution and for f{vx, t) are shown in Fig. [51 (the 
other velocity components show the same behavior). We 
observe three main relaxation stages: 

1. for t < 0.5, a rapid heating is observed which is due 
to acceleration and build up of binary correlations 
in the initially random (uncorrelated) particle sys- 
tem. This is typical for any rapid change of the 
interparticle forces, and proceeds on scales of the 
order of the correlation time, e.g. [21, 22, 23]. 

2. for 0.5 < t < 1.3, the kinetic energy increase sat- 
urates and cooling starts. This means, correlation 
build up is finished and dissipation due to neutral 
gas friction dominates the behavior. 

3. for t > 1.3, the mean kinetic energy decreases ap- 
proximately exponentially, i.e. (£'kin)(i) « e^^'''* 
where the decay constant is found to be 7 « 0.65 « 

The behavior on the third stage resembles a single (Brow- 
nian) particle in a dissipative medium where 7 is the 
velocity relaxation rate corresponding to a relaxation 
time of trei — — 1.54. In case of Brownian parti- 
cles, the velocity distribution rapidly relaxes towards a 
Maxwellian for t < trei- The velocity distributions for 
the present system are shown for four different times 
in Fig. m parts a)-d). The solid curves indicate the 
best fit to a Maxwellian, the obtained "temperatures" 
are shown in Fig. [8] e) by the crosses. The evolution 
towards a Maxwellian is evident which is established 
around t = 2.5. 

This allows us to conclude that, after an initial stage 
(phases 1 and 2), the cluster has reached an equilibrium 
velocity distribution and the subsequent cooling process 
ultimately leading to freezing into a spherical Yukawa 
crystal is well described by local thermodynamic equilib- 
rium: the time-dependent velocity distribution is given 
by fiv,t) ^ exp{-5j^} with kBTit) = 2{Ek^n){t) /i- 
Thus, the system evolves from one equilibrium state to 
another which differ only by temperature. 
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FIG. 8: a) - d) Velocity distribution function f{vx,t) for dif- 
ferent times [as indicated in e)] for A'^ = 27, k = 0.6, v = 3.2. 
Solid lines show the best Maxwellian fit. e) Averaged kinetic 
energy as a function of time. Crosses denote averaged kinetic 
energy obtained from best fit using the equipartition theorem. 



D. Application of Equilibrium Theories to the 
probabihty of metastable states of Yukawa balls 



Based on the results of Subsection lV C[ we expect that 
equilibrium methods such as Monte Carlo or our analyt- 
ical model are applicable to the third relaxation stage. 
Thereby one has to use the equilibrium result for the cur- 
rent temperature T{t). Using temperature dependent re- 
sults such as in Fig.[6l allows one to reconstruct the time- 
dependence of various quantities from the known dynam- 
ics of the kinetic energy: T{t) = r(t^ei)e"^^^*"*"'^- 

Now, the key point is that this local (time-dependent) 
Maxwellian is established long before the particles start 
to feel the potential energy U of the trap and of the pair 
interaction. For example, at t « treii the temperature 
is around 0.15 which is about a factor 100 higher than 
room temperature and one order of magnitude higher 
than the freezing point. In case of very rapid cooling be- 
yond the freezing point the particles will settle (with a 
certain probability) in the stationary state "s" and will 
not have time to escape it since further cooling removes 
the necessary kinetic energy (i.e. the escape probability 
will be low). This means that the decision about what 



TABLE VI: Comparison of experimental results for N = 27 
and A'' = 31 with MD and MC simulations (MC results are for 
the two temperatures T = 0.02 and T = 0.04). Also shown 
are the results of the analytical model ("AM") for T = 0.02 
and with the Boltzmann factor being neglected (T oo). 
For TV = 27 {N = 31) the simulation results are shown for 
K = 0.6 (k = 0.8). 



iV = 27 


P(24, 3) 


P(23, 4) 


P(25, 2) 


Exp sriiiiGiit 


0.46 ± 0.14 


0.46 ± 0.14 


08 ± 06 


MD 


0.46 


0.53 


0.01 


MC(0.02) 


0.56 


0.43 


0.01 


MC(0.04) 


0.43 


0.45 


0.04 


AM(0.02) 


0.67 


0.33 


0.00 


AM(c5o) 


0.12 


0.64 


0.24 


iV = 31 


P(27, 4) 


P(26,5) 


P(25,6) 


Experiment 


0.35 ±0.10 


0.62 ±0.13 


0.03 ± 0.03 


MD 


0.30 


0.59 


0.11 


MC(0.02) 


0.40 


0.55 


0.04 


MC(0.04) 


0.33 


0.50 


0.14 


AM(0.02) 


0.44 


0.53 


0.03 


AM(oo) 


0.02 


0.60 


0.38 



stationary state the system will reach is made at a time 
when the system temperature is close to the melting tem- 
perature. 

Using this idea we compute the probability of 
metastable states from Monte Carlo for two temperatures 
T = 0.02 and T = 0.04, cf. Fig. [6] (at the higher temper- 
ature, due to intershell transitions, shell configurations 
can be identified only with an error of about 8%). We 
also calculate the probability at T = 0.02 within the ana- 
lytical model. Finally we consider the high-temperature 
limit which is obtained by neglecting, in the probabil- 
ity ratios, the Boltzmann factor. The corresponding re- 
sults are presented in Table IVII The overall agreement 
with the experiment is much better than the results for 
room temperature which confirms the correctness of the 
above arguments. Evidently, the Boltzmann factor is cru- 
cial and cannot be neglected, cf. last lines in Table IVll 
The best results are observed for temperatures around 
T — 0.04 which is about two to three times higher than 
the melting temperature where the system is in the mod- 
erately coupled liquid state. This shows that it is indeed 
possible to predict, at least qualitatively, the probabil- 
ities of metastable states in dissipative nonequilibrium 
Yukawa crystals within equilibrium models and simula- 
tions. This is possible in the overdamped limit as is the 
case in dusty plasmas. 



VI. DISCUSSION 

In summary we have presented simulation results for 
Yukawa balls with three different numbers of particles 
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and a broad range of screening parameters and damp- 
ing coefficients. It was shown by extensive molecular 
dynamics and Langevin dynamics simulations that the 
cooling speed (damping coefficient) strongly affects the 
occurrence probabilities of metastable states even if the 
interaction and the confinement remain the same. This 
is similar to the liquid solid transition in macroscopic 
systems where rapid cooling may give rise to a glass-like 
disordered solid rather than a crystal with lower total en- 
ergy. The same scenario is also observed in the present 
finite crystals. While slow cooling leads predominantly to 
the lowest energy state, strong damping gives rise to an 
increased probability of metastable states. These states 
may have an up to five times higher probability than the 
ground state, which is fully consistent with the recent 
observation of metastable states in dusty plasma experi- 
ments [ll|. These metastable states are not an artefact 
of an imperfect experiment or due to fluctuations of ex- 
perimental parameters, but are an intrinsic property of 
finite Yukawa balls. 

Furthermore we showed that screening strongly alters 
the results compared to Coulomb interaction. Generally 
increased screening leads to a higher probability of states 
with more particles on inner shells due to the shorter 
interaction range. An analytical theory for the ground 
state density profile of a confined one-component Yukawa 
plasma also showed that decreasing the screening 

length (increasing k) leads to a higher particle density in 
the center of the trap, which would correspond to a higher 
population of inner shells in our case. 

We presented an analytical model based on the canon- 
ical partition function and the harmonic approximation 
for the total potential energy. This model allowed for 
a physically intuitive explanation of the observed high 
probabilities of metastable configurations. The Boltz- 
mann factor (which always favors the ground state rel- 
ative to higher lying states), competes with two factors 
that favor metastable states: the degeneracy factor [fa- 



voring states with more particles on the inner shell(s)] 
and the local curvature of the potential minimum. Low 
curvature (low eigenfrequency) corresponds to a broad 
minimum and a large phase space volume attracting par- 
ticles. Among all normal modes the dominant effect 
is due to the energetically lowest modes. The thermo- 
dynamic results from Monte-Carlo simulations and the 
analytical theory are in reasonable agreement with each 
other, at low temperatures, as expected. For higher tem- 
peratures anharmonic effects such as barrier heights will 
be equally important. 

It was shown that in thermodynamic equilibrium the 
abundances of metastables are much lower than observed 
in the dusty plasma experiments at the same tempera- 
ture. The reason is that, in equilibrium, the particles are 
given infinitely long time to escape a local potential min- 
imum and they always will visit the ground state more 
frequently than any metastable state. In contrast, in the 
limit of strong damping the particles are being trapped 
in the first minimum they visit. Thus the decision about 
the final stationary state is being made early during the 
cooling process, when the temperature is of the order of 
two to three times the melting temperature. Therefore, 
equilibrium theories without dissipation may be success- 
fully applied to strongly correlated and strongly damped 
nonequilibrium systems. A systematic derivation from 
a time-dependent theory is still lacking and will be pre- 
sented in a forthcoming paper. 
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